function [ output subjectinfo ] = createErsp( subjectinfo, cfg )
%CREATEERSP Create the ERSP for the given subject and sets of
%channels/frequencies.
%
%Usage: Create a cfg with cfg=createErsp_defaultCfg(), edit the cfg and run
%the function with it.
%Output: Contains a report structure that describes this task. (TODO..)


%% Ensure valid output directory
if ispref('ECoG', 'Results')
    resultsBasepath = getpref('ECoG', 'Results');
    basepath = [resultsBasepath 'ERSP\'];
    ensurePath(resultsBasepath);
else
    error('Preference does not contain results path, edit and run SetDefaultPreferences.m.');
end

%% Record the timestamp of the .m file
mfnp = dir([mfilename('fullpath') '.m']);
disp(['Running: ' mfnp.name ' (version ' mfnp.date ')']);
output.mfile = mfnp;

%% For every epoch dataset
eNames = fieldnames(subjectinfo.EpochData);
eSets = {};
for n = eNames'
    nn = cell2mat(n);
    if ~isempty(cfg.epochSets) && ~(any(strcmp(cfg.epochSets, subjectinfo.EpochData.(nn).name)) || any(strcmp(cfg.epochSets, nn)))
        continue;
    end
    eSets{end+1} = subjectinfo.EpochData.(nn);  %#ok<AGROW>
end
if isempty(eSets)
    if isempty(cfg.epochSets)
        disp('There are no epoch sets available for this subject, or they aren''t registered in the SubjectInfo object.');
    else
        warning('The requested epoch sets don''t exist or aren''t available.');
        return;
    end
end
for eSet = eSets
    epochSet = cell2mat(eSet);
    disp(['==============> Creating ERSP''s from epochs: ' epochSet.name]);
    disp(epochSet);
    
    % load the data..
    EEG = pop_loadset('filename', epochSet.dataFilename);
    EEG = eeg_checkset( EEG );

    
    %% run a plot for all requested channels
    % Find channel indices for requested channels
    X = {EEG.chanlocs(:).labels};
    if isempty(cfg.channels)
        IdxTF = true(size(X));
    else
        IdxTF = false;
        for c = cfg.channels
            TF = strcmp(c, X);
            if any(TF)
                IdxTF = IdxTF | TF;
            else
                warning(['Channel "' cell2mat(c) '" was not found!!']);
            end
        end
    end
    channelIdx = find(IdxTF);
    % Enumerate over all indices
    for i = channelIdx
        % Draw the plot
        fig = figure();
        [ersp,itc,powbase,times,freqs,erspboot,itcboot] = ...
            pop_newtimef( EEG, 1, i, [min(cfg.timelimits) max(cfg.timelimits)], [3 0.5] , ...
            'freqs', cfg.frequencies, ...
            'plotitc', 'off', ...
            'plotphase', 'off', ...
            'baseline', [min(cfg.timelimits) max(cfg.timelimits)], ...
            'basenorm', 'on', ... % normalize baseline in the power spectral average instead of dividing by the average power across trials at each frequency (gain model)
            'ntimesout', (max(cfg.timelimits)-min(cfg.timelimits))/10, ...
            'title', [subjectinfo.Name ' - ERSP over ' int2str(EEG.trials) ' trials of ' epochSet.name ' on electrode ' X{i}], ...
            'padratio', 1, ...
            'erspmax', 4); % Set the ERSP clim to [-4 4]
        
        % Save the plot in various formats
        frequenciesStr = [num2str(min(freqs)) '-' num2str(max(freqs)) 'Hz'];
        saveFigure(fig, basepath, [subjectinfo.PrimaryPrefix '-' X{i} ' ' epochSet.name ' ' frequenciesStr]);
        save([basepath subjectinfo.PrimaryPrefix '-' X{i} ' ' epochSet.name ' ' frequenciesStr '.mat'], 'ersp', 'itc', 'powbase', 'times', 'freqs', 'erspboot', 'itcboot');
        close(fig);
    end

end

%% store information in subjectinfo
subjectinfo.Save();
end